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ar-synuclein (aS) is a natively unfolded pre-synaptic protein found in all Parkinson's disease patients as the 
major component of fibrillar plaques. Metal ions, and especially Cu(II), have been demonstrated to 
accelerate aggregation of aS into fibrillar plaques, the precursors to Lewy bodies. In this work, copper 
binding to aS is investigated by a combination of quantum and molecular mechanics simulations. Starting 
from the experimentally observed attachment site, several optimized structures of Cu-binding geometries 
are examined. The most energetically favorable attachment results in significant allosteric changes, making 
aS more susceptible to misfolding. Indeed, an inverse kinematics investigation of the configuration space 
uncovers a dynamically stable /'-sheet conformation of Cu-aS that serves as a nucleation point for a second 
/'-strand. Based on these findings, we propose an atomistic mechanism of copper-induced misfolding of aS 
as an initial event in the formation of Lewy bodies and thus in PD pathogenesis. 
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Parkinson's disease (PD) is the most prevalent neurodegenerative movement disorder affecting more than 
0.1% of the population older than 40 years of age 1 . At the cellular level, PD is characterized by the accretion 
of Lewy body plaques in neurons of the substantia nigra 2 . Lewy bodies are accumulations of amyloid-like 
fibrils, which are comprised primarily of an oligomerized form of a-synuclein (aS) 3 , a natively unfolded protein 
consisting of 140 amino acids 4 . In addition to being implicated in PD, aS also plays a key role in other neuro- 
degenerative disorders, such as dementia associated with Lewy body disease, diffuse Lewy body disease and 
multiple system atrophy 4 . Furthermore, Lewy bodies are found in a significant fraction of Alzheimer's disease 
(AD) patients 5 , suggesting that aS is also important in the pathogenesis of AD. 

An increasing amount of evidence now points to the significance of metal ions in PD. Elevated levels of iron and 
zinc were found in PD brains 6 and high levels of copper were detected in the cerebrospinal fluid 7 . Furthermore, aS 
binds metal ions which accelerate its oligomerization into fibrils 6 . The link between metals and PD is similar to 
that in AD, where altered levels of metals in the brain are also observed 8 . For the case of AD, Bush and Tanzi 9 have 
postulated the metal hypothesis, which argues that neurotoxic effects of proteins are promoted and/or caused by 
their interactions with metal ions. Based on this hypothesis, metal chelation therapies for AD were developed and 
are presently being tested in clinical trials 10 . Similarities between the roles of metals in AD and PD suggest that 
metals represent a promising therapeutic target for PD. 

Copper is one of the most studied metals in neurodegenerative diseases. It appears to have distinctly different 
biological roles in different diseases: while it is harmful and causes structural damage to proteins in AD and PD 11 , 
it is thought to have neuroprotective effects in prion diseases 121415 . In fact, our previous investigation 12 of copper 
attachment to prion protein (PrP) showed that copper binding is a part of PrP's normal function and that this 
protein performs a buffering role, thereby protecting proteins such as aS from damaging effects of free copper 
ions. However, it appears to have detrimental effects in PD, since it was found to be the most effective metal in 
causing the in vitro oligomerization of aS 16 . 

In this report, we use a combination of quantum and classical simulations to study Cu(II) attachment to 
aS. Quantum simulations are employed to determine the local binding geometry of the Cu-aS complex, while 
classical molecular dynamics is applied to probe the structure of the full-length copper-bound aS. These calcula- 
tions reveal that copper attachment to aS causes development of secondary structure in the otherwise unstruc- 
tured protein, which accelerates the misfolding of aS. Based on this finding we propose an atomistic mechanism of 
copper-induced misfolding of aS. 
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Figure 1 | Relaxed 3N10 binding geometry of copper to cr-synuclein. 

The bond distances are: Cu-N s 1.95 A, Cu-Nl 2.01 A, Cu-N2 1.93 A, and 
Cu-O 2.08 A. 



Results 

We consider the attachment of Cu at the experimentally observed 
site His- 50, identified by Rasia et al. 17 using a combination of electron 
paramagnetic resonance (EPR) and nuclear magnetic resonance 
(NMR) techniques. Their data showed that the neighboring residues 
(49-52) are also involved in copper binding. A contribution of the 
N-terminal residues (3-9) was also detected, but a more recent 
study 18 found that the N-terminal residues form an independent 
binding site not related to the primary binding at His- 50. While 
the NMR data indicate good correlation with Cu(II) attachment to 
two nitrogen and two oxygen atoms (2N20), other binding types are 
not precluded. In particular, a 3N10 type of attachment is plausible, 
due to the similarity of its nuclear spin interaction coefficients to 
those of 2N20 19 . Therefore, we consider both 2N20 and 3N10 types 
of attachment in this work. 

To investigate the properties of the copper-aS complex, we con- 
struct multiple binding models of 2N20 and 3N10 attachments. The 
aS is represented by residues 49-54 (VHGVAT) and the copper ion is 
assumed to be coordinated by N^ in the histidine imidazole and by a 
combination of deprotonated amide nitrogens and carboxyl oxygens 
in the protein backbone. Several combinations exist for each binding 
mode, but many can be ruled out because they distort the fragment 
too severely. Using insight from our previous investigations of Cu- 
bound proteins 12 ' 20 ' 21 and results from Ref. 13 , where many possible 
copper attachments modes were investigated, we narrow our search 
to three plausible candidates: Two 2N20 structures, denoted as 
NNOO and NONO according to the order of copper ligands in 
counterclockwise direction starting from the imidazole N^ (see 
Fig. 2), and one 3N10 structure. In the calculations, we use neutral 
ends, i.e., the protein backbone is continued in both directions so that 
it ends with neutral CH3 groups. Each binding model is geometric- 
ally optimized using the hybrid Kohn-Sham/orbital-free density 



Table 1 


Relative energies of copper bine 


ing modes in aS 


Structure 


AE (eV) 


AE (kJ/mol) 


3N10 


0.00 


0.00 


NNOO 


+ 1.77 


170.78 


NONO 


+2.51 


242.18 



functional theory (KS/OF DFT) method 20 and the relaxed geometries 
are shown in Figs. 1 and 2. For all examined structures, the binding 
is mostly planar, with the imidazole ring slightly tilted from the 
coordination plane. Note that for all the binding modes, the 
attachment of copper causes the protein backbone to deflect at 
the binding site. 

Since the 3N10 and 2N20 binding modes differ in the total 
charges (0 and +1, respectively), we include a hydronium ion in 
the solvation shell of the 3N10 model to facilitate a direct compar- 
ison of energies. The relative energies with respect to the 3N10 mode 
are listed in Table 1, which shows that the 3N10 type of binding is 
energetically favorable over the best of the 2N20 structures by 
1.77 eV. 

The properties of full-length aS are studied using classical mole- 
cular dynamics (MD) via NAMD 26 . The model of free aS, whose 
structure has not been determined yet due to its unfolded character 27 , 
is built and equilibrated as described in the Methods Section. 
Extended molecular dynamics simulations show that the free aS is 
mostly linear and extended. No propensity for forming secondary 
structure is observed with the exception of several intermittent turns. 
A snapshot of the equilibrated free aS is shown in Fig. 3 (top). Note 
that the imidazole group in His- 50, which anchors the copper ion in 
copper-bound aS, is extended away from the backbone and can 
readily bind metal ions. 

The structural model of the full copper-bound aS is based on 
the equilibrated free aS, onto which the quantum-mechanically- 
calculated 3N10 copper binding geometry is imposed. Within 
classical molecular dynamics we represent the copper binding by a 
set of harmonic springs. A detailed description of this model is given 
in the Methods section. The initial structure of copperbound aS is 
shown in Fig. 3 (bottom). Note that copper binding causes curving of 
the aS backbone at the binding site, creating a turn-like structure that 
brings fore and aft portions of the backbone into proximity. These 
regions stay proximal during the entire simulation (25 ns). This 
conformation is further stabilized by contacts between side chains, 
the most stable of which are between Glu-46 and Lys-60. A snapshot 
exhibiting this contact is displayed in Fig. 4 (top). The simulation also 
shows that the proximal parts of the protein start to align, indicating 
the development of a /?-sheet around the binding site and the cre- 
ation of a ^-hairpin there. Unfortunately, due to the large time scales 
involved, the formation of sheets in full-length proteins cannot be 
studied in atomistic simulations at present. Instead, we develop mul- 
tiple plausible candidates for a /?-sheet located in the area where 




Figure 2 | Relaxed 2N20 binding geometries of copper to ar-synuclein: NONO (left) and NNOO (right). The bond distances in counterclockwise order 
are: Cu-N.5 1.99 A, Cu-O 2.02 A, Cu-N 1.92 A, and Cu-O 2.18 A (NONO); Cu-N^ 1.98 A, Cu-N 1.97 A, Cu-O 1.93 A, and Cu-O 2.08 A (NNOO). 
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Figure 3 | Misfolding events in or-synuclein due to copper binding. Top: Snapshot of the free a-synuclein showing the location of His-50 (upper). 
Bottom: Initial relaxed structure of Cu(II) -bound a-synuclein. The inset shows the details of copper binding. 



alignment is observed, using inverse kinematics 28 and testing their 
dynamic stabilities in extended MD simulations. Among all the 
candidate structures, only one turns out to be stable during the 
MD simulations. The other ones quickly relax to the unstructured 
form. The dynamically stable sheet conformation involves con- 
tacts between residues 43-45 and 57-59, see Fig. 4 (bottom). It 
remains stable for the entire duration of the simulation (23 ns). In 
addition to the dynamical stability of this sheet, an additional 
structural development at residues 40-42 is also observed during 
the simulation. These residues align with the already existing /?-sheet, 
indicating the development of a second strand parallel to the ori- 
ginal /?-sheet. The observed onset of the formation of a p-sheet has 
important implications, since sheets are a defining characteristic of 



aS in Lewy bodies 4 and this transformation may constitute the first 
step in the development of PD. 

Discussion 

Our calculations show that 3N10 type of copper attachment is pre- 
ferred over 2N20 by 1.77 eV. This conclusion seemingly contradicts 
the work of Rasia et al. 17 , who proposed that their EPR measurements 
indicate a 2N20 copper attachment. However, their analysis of EPR 
data does not include the 3N10 mode; only 4N, 2N20 and 40 were 
considered. Furthermore, since the values of tensor coefficients A|| 
and g|| for 2N20 and 3N10 modes are interspersed 19 , distinction 
between these modes is difficult. To facilitate comparison with the 
experimental results, we calculate the EPR tensor coefficients for 




Figure 4 | Top: Snapshot of evolved Cu-bound cr-synuclein, showing hydrogen bonds between side chains that stabilize the turn-like conformation. 

Bottom: Snapshot of Cu-bound a-synuclein with dynamically stable /?-sheet. For visual clarity, non-participating side chains are not shown. Insets 
show close-ups at the /?-sheet nucleation point. 
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both 2N20 and 3N10 attachments using ORCA 22 . Unfortunately, 
the ORCA EPR calculations do not have sufficient accuracy to dis- 
tinguish between configurations with similar EPR tensors, because 
the calculated tensors typically differ by about 10% from the experi- 
mental results 25 . Our ORCA calculations obtain similar results for 
the2N20 and 3N10 configurations: |A|||=473.8 MHz, g N =2.20 for 
2N20 (NNOO) and |A||| =467.7 MHz and g\\=2.17 for the 3N10 
structure. These values differ from the experimental data of Rasia 
et al. 17 , |A||| = 557.6 MHz and g|| = 2.22, by about 15 %, which is 
similar to the differences between theory and experiments found in 
other works 25 . 

The classical molecular dynamics investigation of the full copper- 
bound aS indicates that the copper attachment in the 3N10 mode 
leads to the development of sheets in aS by bringing the portions of 
the protein around the binding site into close proximity. Due to the 
high computational cost, we did not investigate the 2N20 attach- 
ments. However, since the NNOO mode, which is preferable over the 
NONO structure, has a backbone deflection very similar to that of 
3N10, it will likely follow a very similar structural development path. 
On the other hand, the NONO mode induces a much smaller back- 
bone deflection in the full protein and thus its structural develop- 
ment, if any, will likely be different. 

Based on findings in this work, we propose a mechanism for 
copper-induced misfolding of aS: In the first step, copper attachment 
causes deflection at the binding site and creates a stable turn-like 
structure. This feature is crucial for further structural development, 
since turns play an essential role in the initial steps of protein fold- 
ing 29 . Its presence creates favorable conditions for /?-sheet nucleation 
as adjacent parts of the protein are brought into proximity. Hydrogen 
bond contacts between the proximal side chains are formed and 
further stabilize this conformation. Subsequently, a sheet develops 
around the binding site, serving as a template for the formation of 
additional ^-sheets. This process continues until the protein is fully 
misfolded. 

In summary, we performed computer simulations to investigate 
interactions between a copper ion and aS at an atomistic level. Our 
work, which builds on recent experimental investigations 17 , con- 
sidered 2N20 and 3N10 types of attachment and found that the 
latter is energetically preferred. Further investigation of the preferred 
3N10 mode shows that the copper attachment induces deflection of 
the protein backbone at the binding site, leading to the creation of a 
turn -like feature and backbone alignment as a precursor in the 
development of a ft- sheet. This indicates that the role of copper in 
the fibrillation of aS is both to initiate misfolding via conformational 
modification and to stabilize the partially folded intermediate once 
formed. After the formation of the initial /?-sheet, the folding process 
continues by aligning adjacent parts of the protein to the existing 
/?- sheet, leading to a full misfolding. 

Methods 

The copper binding geometries are investigated at the ab initio level using a hybrid 
DFT/DFT method 20 . In this approach, the chemically active parts of the system, 
including the first solvation shells, are treated at a Kohn-Sham DFT level, while an 
approximate orbital-free DFT is used for the remaining solvent molecules. The 
compatibility between the KS and OF DFT methods enables seamless integration 
between the two. In particular, the flow of solvent molecules across the KS/OF 
interface is allowed and the total energy is conserved. The implementation of the 
hybrid KS/OF DFT approach is based on the RMG 30 " 32 code, in which Kohn-Sham 
equations are solved in real- space and the multigrid technique is used to accelerate 
convergence of the ground state wavefunctions. This method was already successfully 
used to study the binding of copper to the prion protein 12,20,21 . 

In the KS/OF DFT simulations, aS is represented as a fragment consisting of 
residues 49-54 (VHGVAT). Neutral ends are applied to the fragment, resulting in a 
total of 9 1 atoms. In each calculation, the protein fragment is fully solvated and 
minimized using classical molecular mechanics prior to performing KS/OF DFT 
simulations. This is followed by minimization within the KS/OF DFT method, in 
which the protein fragment, the copper ion and 12 water molecules closest to the 
binding site are treated by KS DFT, while 1788 waters are enclosed in the orbital-free 
DFT region. A grid with spacing of 0.32 Bohr is used to represent the wavefunctions, 



corresponding to a kinetic energy cutoff of 48 Ry. Ultrasoft pseudopotentials 31,33 and 
generalized gradient approximation in the PBE form 34 are employed. 

An attempt to differentiate conformational binding modes based on EPR data was 
made using the DFT-based code ORCA 22 . Due to a high computational cost, only the 
immediate binding atoms and their nearest neighbors were included. The dangling 
bonds were passivated by hydrogen. Orbitals were computed via the unrestricted 
Kohn-Sham (UKS) method using the B3LYP 23 exchange-correlation functional with 
the Ahlrichs-TZV basis set 24 , which produced the best results in previous investiga- 
tion of EPR coefficients of copper complexes 25 . The EPR tensor was calculated for the 
copper and contributions from the isotropic, dipole and quadrupole hyperfine 
interactions were all included. 

The simulations of the full-length aS protein are performed with the molecular 
dynamics package NAMD 26 with CHARMM27 35 force-field parameters. As is well 
known, copper and other transition metals are difficult to describe with force fields, 
due to the importance of quantum effects. In this work, the effects of copper were 
modeled in two ways: In the initial model, the copper ion was not present and its effect 
is modeled by fixing its ligands, while we later developed a more sophisticated model 
in which the copper was represented by a set of springs which kept the copper 
geometry close to that of our quantum investigation. Four springs were used to 
represent copper-ligands bonds and two additional springs were used to enforce 
N<5-Cu-N 2 and O-Cu-Ni angles, where the atom labels are from Fig. 1. The angular 
springs are necessary to preserve the planar character of the copper attachment. A 
value of 250 was used as a spring constant for all bond and angular springs. The spring 
constants' units are kcal/(mol A 2 ) and kcal/(mol A) for bond and angular springs, 
respectively. This value reproduced the vibrational amplitudes of the quantum cal- 
culations. The springs were implemented using the tclforces framework available in 
NAMD, and the calculated forces were added to those evaluated from the CHARMM 
force field. Since the overall charge of the 3N10 copper attachment is 0, no charge was 
assigned to the copper ion. The differences between the jS- sheet formation processes 
obtained from the models were very minor. 

The calculations use periodic boundary conditions and solvent is added to the 
simulation cell so that the distance between any protein atom and an edge of the water 
box is at least 7 A. This results in a total atom count that varies between 275,000 and 
350,000. The solvent is ionized with 0.15 mM/L Na + and Cl~ to approximate the 
human intraneuronal ionic concentration 36 and normalize the system to zero net 
charge. The temperature is maintained at 3 10 K by means of Langevin dynamics. The 
Nose-Hoover Langevin Piston algorithm 37 is used to keep the pressure constant at 
1 atm. The Verlet-I/r-RESPA 38 method is employed with the following time steps: 4 fs 
for the long-range electrostatic forces, 2 fs for short-range nonbonded forces, and 2 fs 
for bonded forces. Bonds involving hydrogen are constrained using SHAKE 39 to 
enable the use of large time steps. 

The initial structural model of aS is built according to its amino acid sequence as a 
jS- strand using PyMol 40 . The protein is first equilibrated by running molecular 
dynamics during which the initial structure rapidly degenerates to a conformation 
with almost no secondary structure, with only random coils and transient turns 
present. This equilibration is computationally demanding due to the unfolded nature 
of aS, and requires a large number of solvent molecules. During the 21 ns of simu- 
lation time that could be afforded, the root mean square deviation (RMSD) of the 
protein backbone does not completely flatten out, but other quantities, such as the 
kinetic and potential energies, solvent accessible surface and radius of gyration 
achieve stable values, indicating that the structure of the protein is not changing 
significantly and is sufficiently equilibrated for the study of Cu binding. 

The exploration of possible jS-sheet registrations is aided by inverse kinematics as 
implemented within ProteinShop 28 and based on the AMBER 41 force field. By using 
inverse kinematics we are able to reduce the sample set to contain only energetically 
favorable jS-sheet conformations. 
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